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Abstract 

A brief introduction to the Python computing environment is given. By solving the master 
equation encountered in quantum transport, we give an example of how to solve the ODE problems 
in Python. The ODE solvers used are the ZVODE routine in Scipy and the bsimp solver in GSL. 
For the former, the equation can be in its complex-valued form, while for the latter, it has to 
be rewritten to a real-valued form. The focus is on the detailed workflow of the implementation 
process, rather than on the syntax of the python language, with the hope to help readers simulate 
their own models in Python. 



I. INTRODUCTION 



Python is a general purpose, high-level programming language. It is well documented 
and easy to learn. With its concise and high-readable code, it greatly improves the efficiency 
of code development and code reuse, so we can save our time by simulating our model in 
Python. 

In computation. Python provides us with an easy-to-use environment which saves both the 
code development time and the code execution time. Though it is an interpreted language, 
which means its execution speed is lower compared with the compiled language like C and 
Fortran, its performance has been greatly improved by a set of specific packages designed 
for numerical computation. This is done by gluing well behaved numerical libraries written 
in C/C++ and Fortran and allowing users to use them in Python code, thus we can write 
code efficiently in Python and get an execution speed like C and Fortran. Those packages 
includes Numpy [Il[2], Scipy [HIS], Pygsl [3], Cython [3] and various others packages. Numpy 
is based on the compiled LAPACK libraries that is standard for linear algebra computations. 
It provides the data structure 'array' and fast operations on arrays, such as linear algebra, 
Fourier transform and random number generation, which make it easier and faster to handle 
matrix related problems. As the computation is essentially executed in LAPACK, it can run 
almost as fast as in the C code. Based on Numpy, Scipy provides many modules to perform 
the common tasks in science and industry, such as FFT, sparse matrix, statistics, signal 
processing and ODE solvers. The functionality of Numpy and Scipy is similar to Matlab, 
but they are developed to make scientific computing a natural part of Python, rather than 
be a copy of Matlab. Pygsl is a python interface to the GNU Scientific Library (GSL) [5], 
which is an open source C library for numerical computations in science. The models can 
be expresed in Numpy arrays and then by Pygsl we can use the functions in GSL directly as 
if they are Python functions, which greatly simplifies the usage of GSL functions. Cython 
is based on Python, but it allows for static C type declarations and direct calling of C or 
C++ functions, which combines the high productivity of Python with the execution speed 
of C [6J . Cython is used in the development of the powerful computing software sage [7] (or 
sagemath). Besides computation, there are many powerful visualization tools in Python, 
such as Visual [8] for animation and Matplotlib [9] for Matlab-like plot. Other plotting 
tools, such as Mayavi [10] for 3D visulization and gnuplot [11] for scientific plot, can also be 
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used in Python directly. 

Even though there is very few cases where no package exists for performance improvement, 
we still recommend Python as it could save lots of our time in code development, while using 
the compiled languages can take lots of human time and would be exhausting. On the other 
hand. Python allows us to write the time consuming part of the algorithm, usually long loops, 
directly in compiled languages like C/C++ and Fortran, so we can get an optimal balance 
between human time and machine time. This can be done by lots of powerful wrappers, 
such as F2PY [12], which can automatically wrap Fortran codes and make it callable from 
Python, and Weave [13], which makes it possible to written C and C++ codes directly in 
Python codes. 

In this paper, we will use Python to solve the master equation in reference [11], which is 
an ordinary differential equation (ODE). This is an introduction to the Python computing 
environment and a detailed example of solving the ODE problems in Python. We will use 
three solvers, the ZVODE routine in Scipy, the bsimp solver and the rkf45 solver of GSL to 
solve it. We will also give an example of using Cython to optimize the code to get faster 
execution speed. The emphasis is layed on the workflow of how to solve the ODE problems 
in Python rather than on the syntax of this language. One can turn to the books [TSHTT] 
for a detailed tutorial of computational physics in Python. 



II. THE MASTER EQUATION 



The master equation to be solved is 



f)= -i [n, p] - r [p - diag{p)] , 
/ n n n \ 



(1) 



where the Hamiltonian % 



V 



























23 



describes the coherent tunneling via 



adiabatic passage (CTAP) scheme, F is the T2 dephasing rate and VL12 and 1^23 are the 
tunneling rate between the corresponding quantum dots 



{t) = or"'-" exp 

1^23 [t) = Vr"'"" exp 
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This is an evolution equation of the density matrix of the system, which describes the evo- 
lution of the electron in the dots. The master equation is often used when the decoherence 
effect are considered due to the interaction between the system and its environment. The 
diagonal terms of the density matrix describe the population of the electron in the corre- 
sponding dots, while the off-diagonal terms describe the correlation between corresponding 
dots. For the T2 dephasing term, the off-diagonal terms describe the decoherence of the 
electron, which corresponds to the loose of its quantum nature. The neglected diagonal 
terms in the dephasing terms describe the loss of particles in corresponding dots due to the 
interaction with the environment, which is not important in the CTAP scheme. 

To judge the efficiency of the transporting system, this equation needs to be solved. We 
would first solve the case without the dephasing term using the solver ZVODE and bsimp, 
and only consider the case of counter-intuitive pulses, that is, the result shown in Figure. 3(b) 
in reference |T3j. After that, we would solve the case with the dephasing term included using 
the rkf45 solver and give the pseudo-color plot correspoding Figure. 4 in reference [H]. The 
parameters are chosen as the same as there with a = tmax/8. Reducing t by -k/Q"^"'^, that is, 
t = t/ (7r/i7™"^) and tm = tmax/ (vr/i7'"'^'^), and then rescale t to 1 by t = t/tm, the equation 
is reduced to the form 



This is the equation we need to solve in Python. 

III. THE SOLVERS 

Scipy provides two ODE solvers 'odeint' and 'ode' in its 'integrate' module, with the 
former using the 'Isoda' of the Fortran library odepack and the latter using the VODE ( 
for real-valued equations) and the ZVODE (for complex-valued equations) routines. The 
ZVODE routine provides the implicit Adams method for non-stiff problems and a method 
based on the backward odifferentiation formulas (BDF) for stiff problems. Stiff equation is 



P 



mtrn\H,p] - VTxt,n[p - diag{p)], 



(2) 



where F = T/fl"^""^ and the coupling pulses are reduced to 
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that includes some terms that can lead to rapid vibration in the solution and it requires the 
step size to be taken extremely small if using non-stiff solvers. This would happen when 
the variables are changing on two vastly different scales. If you do not know whether you 
problem is stiff or not, using the stiff solvers would be safe, although it would cost more 
time. Here we will use the ZVODE routine with the BDF method. The CTAP scheme is 
non-stiff, but we use stiff solvers here just to give examples of how to use them. In the 
last part when we solve the equation with the dephasing effect, we would use the non-stiff 
solver and optimize it by Cython, as there is a hard requirement on computation time there. 
As the solver only accepts a set of equations that are in a vector form, we have to rewrite 
equation ^ from the matrix form to a one-dimensional array form. This can be done via 
the Numpy array and related operations. The equation finally should be in the form 

^ = -intmf, (4) 

where y is the column version of p and / is the derivative of y with time. 

Listing[l]impliments the solution of equation (|4]). We first import it by 'from Scipy.integrate 
import ode', which imports the solver ode from Scipy's integrate module. Then we define 
equation Q and its Jacobian using the data structure array provided by Numpy, which is 
straightforward in describing matrix related problems . The integrator of the solver should 
be set to use the ZVODE routine with the BDF method. After setting the step size and 
the precision anticipated, we can call this solver to forward the calculation step by step 
until it reaches the final value. We can plot the data directly in Python using Matplotlib 
or gnuplot [18]. Here we output the data to a file and plot it in gnuplot. The source file of 
gnuplot would also be given in listing |2] and the result is shown in Figure. [T] 

Listing 1. This program solves the CTAP regime without the dephasing term using the zvode 
routine of Scipy. 



#! /usr/bin/env python2.7 

#First import the modulus used in this program. Numpy is so often used that 
perhaps every 

#program should import it first. The ode solver is in Scipy's integrate 
modulus . 
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import numpy 

from scipy . integrate import ode 

#Define function dy/dt = f. We have to write the elements of f out by hand, 
as 

#this solver requires the Jacobian , whose elements must be calculated by hand 

#For the rkf45 solver used below, which does not require the Jacobian , we can 
express it by 

#Numpy array operations without calculating its elements by hand. 
def userSupply (t , y, tm) : 

j l=numpy . exp (— 32.*(t— 9./16)**2) 

j 2=iiumpy . exp (— 32.*(t— 7./16)**2) 

#complex number is represented by j in Python. 5j means 5* i and Ij means 
i . 

return — lj*numpy. pi*tm*numpy. array ( 

[jl*y[l]-jl*y[3] , jl*y[0] + j2*y[2]-jl*y[4] , j 2 *y [1] - j 1 *y [ 5 ] , 
jl*y[4]-jl*y[0]-j2*y[6] , j 1 *y [3] + j 2 *y [5] - j 1 *y [1] - j 2 *y [ 7] , 
j2*y[4]-jl*y[2]-j2*y[8] , j 1 *y [7] - j 2 *y [3 ] , 
jl*y[6] + j2*y[8]-j2*y[4] , j 2 *y [7] - j 2 *y [ 5 ] ] ) 

#The jacobian df/dy of the ODE equations . 
def j ac ( t , y , tm) : 

j l=iiumpy . exp (— 32.*(t— 9./16)**2) 

j 2=iiumpy . exp (— 32.*(t— 7./16)**2) 

return — Ij *numpy. pi*tm*numpy. array ([ [0 , jl, 0, — j 1 , 0, 0, 0, 0, 0], 

[jl , 0, j2 , 0, -jl , 0, 0, 0, 0] , 
[0, j2 , 0, 0, 0, -jl , 0, 0, 0] , 
[-jl , 0, 0, 0, jl , 0, -j2 ,0,0], 
[0, -jl , 0, jl , 0, j2 , 0, -j2 , 0] , 
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[0, 0, -jl , 0, j2 , 0, 0, 0, -j2] , 
[0 , 0, 0, -j2 , 0, 0, 0, jl , 0] , 
[0, 0, 0, 0, -j2 , 0, jl , 0, j2] , 
[0, 0, 0, 0, 0, -j2 , 0, j2 , 0]]) 

#initial time and value. 
to = 

yO = numpy . array ([ 1 , 0, 0, 0, 0, 0, 0, 0, 0], dtype=numpy . complexl28 ) 
#t-max and the step size . 
tm = 8000 
H = l.e-3 

#Call the solver ode with the function ' userSupply ' and the Jacobian 'jac ' , 

and 

#define the integrator ' zvode ' with the 'bdf' method and set tolerance 'rtol 

r = ode ( userSupply , jac ). set -integrator ( 

'zvode', niethod=' bdf ' , with_j acobian=True , rtol=l.e — 6, order=5) 

^Prepare the intial values and the parameters t_max . 

r. set_initial_value (yO, tO) . set_f_params(tm) . set_jac_params(tm) 

#define the file to write data to. 
output l=opon ( ' dataCounBdf ' , 'w') 

#Forward the integrator in a loop and output the data to the file. 
while r. successful and r.t < 1.: 

if ( 1. - r.t) < H: H = 1. -r .t 

r . integrate ( r . t + H) 

outputl .write ("%f %f %f %f\n" % (r.t , r.y[0], r.y[4], r.y[8])) 
outputl . close ( ) #close the file. 



Listing 2. The program plots the evolution result, 
set terminal tikz standalone color size 5in,3in 
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solution of the ZVODE routine 
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FIG. 1. Solution got via the ZVODE routine for the counter- intuitive pulse without dephasing. 
The particle evolves from the first dot to the third dot without any population on the second dot. 
This is the result of the CTAP scheme, as explained in reference [Hj. Here tm has been set to 8000 
to avoid vibrations of the lines. 



set out 'counterBdf.tex' 
3 set xlabel ' $\ f r ac { t }{tm} $ ' 

set ylabel '$\rho$ ' 
5 set yrange [—0.1:1.2] 

set grid 

7 set title 'solution of the ZVODE routine' 

plot 'dataCounBdf ' using ($1):($2) smooth csplines title ' $\ rho_ { 1 1} $ ' , ' 
dataCounBdf using ($1):($3) smooth csplines title '$\rho_{22}$ ' , ' 
dataCounBdf' using ($1):($4) smooth csplines title '$\rho_{33}$ ' 

GSL provides many ODE solvers, using the implicit Bulirsch-Stoer method, the Gear 
method and various Runge-Kutta methods. They can be used in Python via the Pygsl 
module, which makes their usage straightforward, just like using the solvers in Scipy, and 
you can read its document by typing 'help('pygsl.odeiv')' in Python. Here the solver bsimp 
is used, which is given in listing [3] The bsimp solver impliments the implicit Bulirsch-Stoer 
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method of Bader and Deuflhard. For smooth functions, the Buhrsch-Stoer method is the 
best way to achieve both high-accuracy solutions and computational efficiency [12], with 
the implicit Bulirsch-Stoer method designed for stiff problems. As above, the Pygsl module 
should to be imported first. Then is the definition of the equation and its Jacobian. It 
should be noticed that the solvers in GSL only accept real-valued equations, so y's real and 
imaginary part in equation Q have to be separated and form a new set of coupled equations 
with 18 elements, with ?/[2i] and ?/[2z + 1] the real and the imaginary part of the original|/[i]. 
Their preparation is a little different from solvers in Scipy, which can be learned from its 
documentation. It only needs inputs of the initial step size and the errors tolerated, as the 
following step sizes are adjusted automatically to optimize its speed and precision. The 
gnuplot source code is as the code [2] except with a different input data name and output file 
name, so we do not put it here. The result is given in Figure. [2] 

Listing 3. This program solves the CTAP regime without the dephasing using the bsimp solver of 
GSL. 
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/usr/bin/env python2.7 
#First import the modulus used in this program. 
#The solver is in pygsl 's odeiv module. 
import nunipy 
from pygsl import odeiv 

#User supplied function dy/dt = f. This is more complex compared with the 
above code . 

#As we have to seperate the real and imaginary parts from the original 

equations . 
def userSupply ( t , y, tm) : 

j2 = numpy. exp( -32.*(t -9./16) **2) 
jl = numpy . exp ( — 32 . * ( t — 7. / 16) * *2) 
return numpy. pi *tm*numpy . array ( 

[jl*(y[3]-y[7]) , j 1 * ( y [6] - y [ 2 ] ) , j 1 * (y [1] -y [ 9 ] ) +j 2 *y [ 5 ] , 
jl*(y[8]-y[0])-j2*y[4] , j 2 *y [3] - j 1 *y [ 1 1 ] , j 1 *y [10] - j 2 *y [ 2 ] , 
jl*(y[9]-y[l])-j2*y[13] , j 1 * ( y [0] - y [ 8 ] ) +j 2 *y [ 1 2] , 



jl*(y[7]-y[3])+j2*(y[ll]-y[15]) , j 1 * (y [2] -y [ 6] ) +j 2 * (y [14] - y [ 1 0] ) , 
j2*(y[9]-y[17])-jl*y[5] , j 2 * (y [16] - y [ 8 ] ) +j 1 *y [4] , 
jl*y[15]-j2*y[7] , j 2 *y [6] - j 1 *y [ 1 4] , j 1 *y [13] + j 2 * (y [1 7] - y [ 9 ] ) , 
-jl*y[12] + j2*(y[8]-y[16]) , j2 *(y[15] -y [1 1] ) , j 2 * (y [10] - y [ 1 4] ) ] ) 

#the Jacobian df/dy and df/dt. The bsimp solver requires both df/dy and df/dt 

def j ac ( t , y , tm) : 

j2 = numpy . exp ( — 32.* (t — 9./16) **2) 
jl = numpy . exp ( — 32.* (t — 7./16) **2) 
dj2dt = -64*(t -9./16)*j2 
djldt = -64*(t -7./16)*jl 
dfdy = numpy. zeros ( [18 , 18] ) 



dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 
dfdy 



0, 3] = jl; dfdy[0, 7] = -j 1 ; 

1, 2] = -jl; dfdy[l, 6] = jl; 

2, 1] = jl; dfdy [2, 5] = j2; dfdy [2 , 9] = -j 1 

3, 0] = -jl; dfdy [3, 4] = -j2; dfdy [3 , 8] = jl 

4, 3] = j2; dfdy[4, 11] = -j 1 

5, 2] = -j2; dfdy[5,10] = jl 

6, 1] = -jl; dfdy [6, 9] = jl; dfdy[6, 13] = -j2 

7, 0] =jl; dfdy [7, 8] = -j 1 ; dfdy [7 , 12] =j2 

8,3]= -jl ; dfdy[8 , 7] = jl ; dfdy [8 , 11] = j2 ; dfdy[8 , 15] = 

9, 2] =jl; dfdy [9, 6] = -j 1 ; dfdy[9, 10] =-j2; dfdy [9, 14] 

10, 5]=-jl; dfdy[10, 9]=j2; dfdy[10, 17]=-j2 

11, 4] = jl; dfdy[ll, 8] = -j2; dfdy[ll, 16] = j2 

12, 7] = -j2; dfdy[12, 15] = jl 
13 , 6] = j2 ; dfdy [13 , 14] = -j 1 

14, 9] =-j2; dfdy[14, 13] = jl; dfdy[14, 17] = j2 

15, 8] =j2; dfdy[15, 12] = -j 1 ; dfdy[15, 16] =-j2 

16, 11] = -j2; dfdy[16, 15] = j2 



-j2 
j2 
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dfdy[17, 10] = j2; dfdy[17, 14] = -j2 
dfdt — numpy . array ( 

[djldt*(y[3]-y[7]) , dj Idt * (y [6] -y [ 2 ] ) , dj Idt * (y [1] -y [ 9 ] ) +dj 2dt *y [5 ] , 
djldt*(y[8]-y[0])-dj2dt*y[4] , dj 2dt *y [3] - dj Idt *y [ 1 1 ] , djldt*y[10]- 
dj2dt*y[2] , 

djldt*(y[9]-y[l])-dj2dt*y[13] , dj Idt * (y [0] - y [ 8 ] ) +dj 2dt *y [ 1 2] , 
djldt*(y[7]-y [3])+dj2dt*(y[ll]-y[15]) , dj 1 dt * (y [2] -y [ 6 ] ) +dj 2dt * (y 

[i4]-y[io]) , 

dj2dt*(y[9]-y[17])-djldt*y[5] , dj 2dt * (y [16] - y [ 8] ) +dj Idt *y [4] , 
djldt*y[15] - dj2dt*y [7] , dj2dt*y[6] - djldt*y [14] , dj Idt *y [13] + dj2dt * (y 

[i7]-y[9]) , 

-djldt*y[12] + dj2dt*(y[8]-y[16]) , dj 2dt * (y [15] - y [ 1 1 ] ) , dj 2dt * (y [10] - y 
[14])]) 

return numpy . pi *tni*dfdy , numpy . pi *tm* dfdt 

#Initial values . 
dimension = 18 
t = 

y = numpy. array ([1 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) 
tm = 50 

#the stepping function advances the solution from t-i to t-i+1 . 
#step-bsimp is the bsimp solver we use. 

step = odeiv . step.bsimp ( dimension , userSupply , jac , args=tm) 
#the control function optimizes the step size with the input tolerated 
errors . 

control = odciv . control_yp_new ( step , 1 . e — 6 , 1 . c — 6) 
#Based on the stepping and the control funtion , the evolve 
#function advances the solution in a given interval . 
evolve = odeiv . evolve ( step , control , dimension ) 
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solution of the solver Bsimp 
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FIG. 2. Solution got via the bsimp solver of GSL for the counter-intuitive pulse without dephasing. 
It's the same as Figure. [TJ tm also equals to 8000. 
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#The file to store data. 

output l==open ( ' datalnBsipmp ' , 'w ' ) 

#initial step size. It can optimize step size under given errors. Only need 

to input an initial step size, 
h = l.e-2 
while t < 1 . : 

if (1. - t) < h: h = l. - t 

t, h, y = evolve . apply (t , 1, h, y) 

outputl . write ("%f %f %f %f\n" % (t, numpy . sqrt (y [0] * *2 + y [ 1 ] * * 2 ) , numpy . 
sqrt(y[8]**2 + y[9]**2) , numpy. sqrt(y[16]**2 + y[17]**2))) 

#clo s e file. 
outputl . close 
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Now we add the dephasing term and solve the equation ([2]). This is done by the rkf45 
solver of GSL, which uses the Runge-Kutta-Fehlberg method with adaptive stepsize control 
using the 5th order error estimate. The rkf45 method is the most general method and applies 
to almost all initial value ODE problems. The solution process gets greatly simplified using 
this solver, as it does not require the Jacobian of the system, the calculation of which 
takes most of our time. The Jacobian needs to do partial differentiation on the funtion / 
of equation (|4]), which needs to be done by hand. Without it, we can use Numpy array 
operations to write the equation ^ to the one-dimensional array form ^ and do not need 
to calculate by hand at all, thus greatly saves our time. This can be done by substituting 
P = Pr + Wc into equation (|2| and seperating the derivative of the pR and pc- After 
rescaling, we would get the following equation 

dpR 



dt 
dpc 



-- [H,pc\ - T[pr - diag{pR)] 

-[H, Pr] - T[pc - diagiypc)]. (5) 



dt 

This form can be expressed to the one- dimensional form easily in Numpy, as the listing |4] 
shows, which solves the equation with dephasing included. 

To get the pseudo-color plot of Figure. 4 of reference [I^. we need to solve the master 
equation ([2]) many times with different values of F and tm each time. Here we solve it 
for 100x500 times with 100 different F and 500 different t^. This is done in a loop and 
would cost huge computation time. We would give an Cython code to optimize the Python 
code, which is also an example of how to use Cython. The optimzed code with Cython is 
given in listing [6] in the next section. The gnuplot source code is given in listing |5] and the 
pseudo-color plot is shown in Figure. |3] 

Listing 4. This program solves the CTAP regime with the dephasing included using the solver 
rkf45 of GSL. 



/usr/bin/env python2.7 
#First import the modulus used in this program. 
#The solver rkf4-5 is in pygsl 's odeiv module. 
import numpy 
from pygsl import odeiv 
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#Abstract class used for declaration . 
class Function : 

def evaluate ( self , t , y , li s ) : 
tm , x= 11 s 
return 

#User supplied function dy/dt = f. We can use the Numpy array and related 
op erations to define it , including array reshape and concatenate , which 
simplify s the process. 
class userSupply ( Function ) : 

def evaluate ( self , t , y , 11 s ) : 
tm , x= 1 i s 

j 1 = numpy . exp (— 32.*(t — 9./16) **2) 

j2 = numpy . exp ( —32.* (t —7./16) **2) 

Ifenumpy. array ([[0 , j 1 , ] , [ j 1 , , j 2 ] , [0 , j2 , ] ] ) 

rhoR=y [0 :9] . reshape ((3 ,3) ) 

rhoC=y [9: 18] . reshape ((3,3)) 

dRdt=numpy . dot (H, rhoC )— numpy . dot (rhoC ,H)— x* (rhoR— numpy . diag ( 
rhoR)) 

dCdt=numpy . dot (rhoR ,H)— numpy . dot (H, rhoR)— x* (rhoC— numpy . diag ( 
rhoC ) ) 

return numpy . pi *tm*numpy . concatenate (( dRdt . reshape (( 9 ) ) ,dCdt. 

reshape ((9) ) ) ) 

#compute-dephase impliments the loop within which the equation is solved with 

different parameters each time. 
def compute_dephase ( f , Nl ,N2) : 

outputl=open ( 'dataDcphascl ' , 'w') 

output2=open ( ' dataDephase2 ' , 'w') 

output3=open ( ' dataDephase3 ' , 'w' ) 

dimension=18 
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for X in numpy . linspace ( . , . 5 ,N1) : 

for tin in numpy . linspace ( , 1 ,N2) : 

step = odeiv . step_rkf45 (dimension , f . evaluate , args=(tm 
, x)) 

control = odeiv . control_yp_new ( step , 1 e — 7,le— 7) 
evolve = odeiv . evolve ( step , control , dimension ) 
y=numpy . zeros ( ( dimension ) ) 
y[0] = 1.0 
t=0.;h = l.e-4 
while t <1.: 

if (l.-t)<h: h=l.-t 
t, h, y = evolve . apply (t , 1 . ,li,y) 
output 1 . write (" %f " % y[0]) 
output2 . write (" %f " % y[4]) 
outputs, write (" %f " % y[8]) 
outputl . write ("\n" ) 
output2 . write ("\n" ) 
outputs . write (" \n" ) 
outputl . close 
output2 . close 
outputs . close 

compute -dephase (userSupply() ,2,S) 

Listing 5. The program plots the pseudo-color plot for the master equation with dephasing. 

set terminal post color enhanced 

set out 'dephasel.ps' 

set xlabel 't_{max}' 

set ylabel '{/Symbol G} ' 

set xtics ("0" 0, "2" 100, "4" 200, "6" S00,"8" 400, "10" 500) 

set ytics ("0" 0, "0.1" 20, "0.2" 40,"0.S" 60, "0.4" 80, "0.5" 100) 
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set cbrange [0:1] 
set pmSd map 

set title '{/Symbol r}_{l}' 
10 splot ' dataDephasel ' matrix 
exit 



IV. CYTHON 

Cython allows C type declarations both for variables and functions, and you can use C 
functions from C libraries in Cython. It can improve the Numpy array operations, especially 
array iteration, to the speed of C array operations, so we can write code efficiently using 
Numpy array while get the speed of C. To use Cython, there are generally three stages. 
Firstly we typedefine the variables to C static types via the key word 'cdef ', which removes 
the dynamical nature of those variables and greatly improves the operation speed involving 
them. If their are Numpy arrays, we should also typedefine their type, which is shown 
in the listing [6j There are a static type and a dynamical type, both of which should be 
typedefined. After this first stage, the program would get a faster speed compared with 
the original one, especially when the iterator in the long loops are typedefined. The second 
stage is to typedefine the Python functions to return static type value via the key word 'cdef 
type-used function (tpye-used variable)', where the type-used are appropriate types. This 
can greatly improve the speed of function calls, but now this function can only be called 
in Cython. If using the key word 'cpdef, then the function can be called in Python. The 
function type defination can improve the behavior of the total code greatly. The last stage is 
to turn off some regular checks for some functions in the code, especially for inline functions. 
Well optimized Cython code can get an speed as fast as C and Fortran, as can be seen in 
the examples provided by the Cython website. For detailed documentation, turn to the 
reference [H E] and references there. The listing [6] below only serves as an example of how 
to use Cython. The list [7] shows the results of profiling the Python code and the Cython 
code. The cpu is '2 Intel(R) Core(TM)2 Duo CPU T7500 @2.20GHz'. We can see that 
the Cython code saves about one third time compared with the Pythond code. The speed 
got is not very much compared with hundreds of times in examples provided by Cython 
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documentation. There are two reasons for this. The first is that this is not optimizing a 
pure Python code. We have used many well behaved packages in the Python code, which is 
very fast already. The second is that our Cython code is not perfect. Many modifications 
can be done to optimize it further. 

Listing 6. The Cythond code to optimize the Python one to get speed up. 
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/usr/bin/env python2 .7 

#First import the modulus used in this program. 
#The solver rkf45 is in pygsl 's odeiv module. 

import numpy as np 
cimport numpy as np 
from pygsl import odeiv 
cimport cytlion 

cdef extern from "matli.li": 

double exp( double) 
@cython . profile (False) 
cdef inline double EXP(double x) : 

return exp(x) 

DTYPE=np. float64 

ctypedef np.float64_t DTYPE.t 

#Abstract class used for declaration . 
cdef class Function: 

cpdef np . ndarray [DTYPE_t , ndim = l] evaluate ( self , double t , np . ndarray [ 
DTYPE_t,ndim = l] y, np . ndarray [DTYPE_t , ndim = l] lis): 
cdef tm,x 
tm , x= 1 i s 

cdef np. ndarray [DTYPE_t,ndim=l] a=np . zeros ( ( 1 8) , dtype=DTYPE) 
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45 



47 
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return a 

#User supplied function dy/dt = f. We can use the np array and related 

operations to define it , including array reshape and concatenate , which 
simplify s the process. 
cdef class userSupply ( Function ) : 

cpdef np . ndarray [DTYPE_t , ndim = l] evaluate ( self , double t , np . ndarray [ 
DTYPE_t,ndim = l] y , np . ndarray [DTYPE_t , ndim = l] lis): 
cdef tni,x 
tm , x= li s 

jl = EXP(-32.*(t-9./16)**2) 
j2 = EXP(-32.*(t-7./16)**2) 

cdef np. ndarray [DTYPE_t,ndim=2] H=np . zeros ( ( 3 , 3 ) , dtype=i)TYPE) 
H[0,l] = jl;H[l,0] = jl;H[l,2] = j2;H[2,l] = j2 
rhoR=y [0 :9] . reshape ((3 ,3) ) 
rhoC=y [9: 18] . reshape ((3,3)) 

dRdt=np . dot (H, rhoC )— np . dot (rhoC ,H)— x* (rhoR— np . diag (rhoR) ) 
dCdt=np . dot (rhoR ,H)— np . dot (H, rhoR)— x* (rhoC— np . diag (rhoC ) ) 
return np . pi*tm*np . concatenate ( (dRdt .reshape((9)) ,dCdt . 
reshape ((9) ) ) ) 

#compute-dephase impliments the loop within which the equation is solved with 

different parameters each time. 
cdef compute.dephase ( Function f , double Nl, double N2) : 

outputl=open ( ' dataDephasel ' , 'w') 

output2=open ( ' dataDephase2 ' , 'w') 

output3=open ( ' dataDophase3 ' , 'w ' ) 

cdef int dimension=18 

cdef double x,tm 

cdef double t , h 



19 



cdef np. ndarray [DTYPE_t,ndini=l] arg=np . zeros ( ( 2 ) , dtype^JTYPE) 

cdef np . ndarray [DTYPE_t , ndim = l] y=np . zeros (( dimension ), dtype=DTYPE) 

for X in np . linspace (0 . ,0 . 5 ,N1) : 

for tm in np . linspace (0 , 10 ,N2) : 
arg [0]=tm; arg [l]=x 

step = odeiv . step_rkf45 (dimension , f . evaluate , args=arg) 
control = odeiv . control_yp_new ( step , 1 e — 7,le — 7) 
evolve = odeiv . evolve ( step , control , dimension ) 
y[0] = 1.0 
t=0.;li = l.e-4 
while t<l.: 

if (l.-t)<h: h=l.-t 
t, h, y = evolve. apply(t ,1. ,li,y) 
outputl .write (" %f " % y[0]) 
output2. write (" %f " % y[4]) 
outputs, write (" %f " % y[8]) 
y=np . zeros ( ( dimension ) , dtype=DTYPE) 
outputl . write ("\n" ) 
output2 . write ("\n" ) 
outputs . write ("\n" ) 
outputl . close 
output2 . close 
outputs . close 

cpdef callFunction ( Function f , double Nl, double N2) : 
compute_depliase ( f ,N1 ,N2) 



Listing 7. The result of profiling the Cython code and the Python one. 



#Profile result for dephaseCython . pyx . 

#The CPU of my laptop is 2 Intel (R) Core(TM)2 Duo CPU T7500 @2.20GHz. 
Tue May S 21:24:15 2011 Profile. prof 



20 



314988757 function calls in 4075.919 seconds 



Ordered by: internal time 



ncalls tottime percall cumtime percall filename : lineno ( function ) 



1 



3 



5 



7 



9 



1 



3 



5 



7 



3987941 


2750 


219 





001 


4032 


149 





001 


{pygsl . ..callback . 




gsl_odeiv 


_e vol ve 


-apply} 












61192534 


693 


108 





000 


1281 


930 





000 


twodim_base . py : 220 ( diag 


) 


61192534 


229 


666 





000 


413 


172 





000 


numeric .py:216( as array) 




61192534 


183 


506 





000 


183 


506 





000 


{numpy. core . multiarray . 


array 


122385068 


175 
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000 
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000 


{Icn} 




1 


23 


031 


23 


031 
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919 


{dephaseCython . callFunction } 


3987941 


17 


043 





000 


4049 


192 





001 


odeiv.py:421( apply) 




50000 





860 





000 


1 


364 





000 


odeiv.py:385( __init__) 




50000 





392 





000 





649 





000 


odeiv.py:88( __init__) 




200000 





318 





000 





318 





000 


{hasattr } 




50000 





287 





000 





476 





000 


odeiv .py:407( ..del..) 




50000 





273 





000 





439 





000 


odeiv.py:260( __del__) 




50000 





263 





000 





447 





000 


odeiv. py:110( ..del..) 




50000 
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000 





317 





000 


odeiv. py:365( ..init..) 




50000 
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000 
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{pygsl . __callback . 
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iv.step 
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{pygsl . ..callback . 
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50000 
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50000 
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{pygsl . ..callback . 





gsl _odeiv_step - free} 
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50000 





110 
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{pygsl . --Callback . 




gsl 


_odc 


iv_cvolve 


_free } 












29 


50000 





091 





000 
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000 


{pygsl . ..callback . 




gs] 


_odeiv_control_free} 














50000 





080 





000 





080 





000 


odeiv.py:163(_get_func) 


31 


50000 





077 





000 





077 





000 


odeiv.py:166( _get_jac) 




50000 





077 





000 





077 





000 


odeiv.py:169( _get_args) 
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50000 





077 





000 
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000 


odeiv.py:160( _get_ptr) 




50000 





076 





000 
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000 


odeiv.py:293( _get_ptr) 
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function_base .py:6( linspace) 
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{method 'disable' of '_lsprof. 






Profiler ' 


obj acts } 













39 



41 



43 



45 



47 



49 



51 



53 



#Profile result for dephase.py. 

Tue May 3 23:33:12 2011 Profile. prof 

651748000 function calls in 6163.858 seconds 

Ordered by: internal time 

ncalls tottime percall cumtime percall filename : lineno ( function ) 



30596267 2708.150 
91788801 1223.705 
61192534 674.102 
122385068 455.155 
122385068 363.905 
ndarray ' objects} 



0.000 5925.743 

0.000 1223.705 

0.000 1259.449 

0.000 455.155 

0.000 363.905 



0.000 dephase . py : 14 ( evaluate ) 

0.000 {numpy. core . multiarray . array} 

0.000 twodim_base.py:220(diag) 

0.000 {numpy. core . multiarray . dot} 

0.000 {method 'reshape' of 'numpy. 
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000 


411 
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000 


numeric . pv :216f as array) 


55 
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{pygsl . __callback . 




gsl- 


odeiv.evolve 


-apply} 
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000 
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000 


{len} 


57 
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000 


99 
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000 


jnumpy. core . multiarray . 




concatenate} 


















1 


16.182 


16 
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dephase .py:26( compute_dephase ) 
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odeiv .py:421( apply) 




50000 
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000 


1 
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000 
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odeiy.py:88( __init__) 
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000 
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000 


{hasattr } 


63 


150300 


0.301 





000 
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000 


{method 'write' of 'file' 



objects } 



50000 


0.278 0.000 
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000 


odeiy.py:407( __del__) 


50000 


0.263 0.000 
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000 


odeiy.py:260( __del__) 


50000 


0.255 0.000 
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odeiv. py:110( __del__) 


50000 


0.214 0.000 
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odeiy .py:365( __init__) 
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{pygsl . ..callback . 
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109 





000 


{pygsl . ..callback . 
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{pygsl . ..callback . 
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_odeiv_step_free} 
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{pygsl . ..callback . 


gsl 


_odeiy_evolve_free} 
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088 





000 


{pygsl . --Callback . 


gsl 


_odeiv_control_free} 
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V. CONCLUSION 

The above examples are about how to solve the ODE problems in Python. One can also 
do monte carlo simulation in Python. There are many packages that can be used in Python, 
such as ALPS (Algorithms and Libraries for Physics Simulations) [201 EI] and PyMC |22]- 
GSL also has modules for Monte Carlo simulation , which can be used via Pygsl as above. 
Python is powerful and easy to learn. Its syntax is simple and most of the time, it is where 
to find the libraries needed and combine them together that provides difficult for beginners. 
Once get familiar with it, it would be elegant and intuitive to do numerical simulations. 
We hope this can help non-computation specialists get familiar with Python and implement 
their own models efficiently in it. 
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